The present invention relates to a method for the simultaneous optimization of an arbitrary number of electromagnetic pulses, which act in a cooperative way, or mutually compensate each other's errors. That is, the invention generally relates to pulses which can have improved properties when cooperating with each other than single pulses. In experiments with several scans, undesired signal contributions can be suppressed by COOP pulses, which complements and generalizes the concept of phase cycling. COOP pulses can be optimized efficiently using an extended version of the optimal-control-based gradient ascent pulse engineering (GRAPE) algorithm. The advantage of the COOP approach is demonstrated experimentally and theoretically for broadband and band-selective excitation and saturation pulses.
It is generally known that radio frequency pulses can influence the system states of spin systems. In this context, according to the present invention, the expression “system state” signifies all physical states a plurality of coupled or uncoupled spins can take. In this context, a system state can be composed of one or several single contributions. Accordingly, Levitt (Prog. Nucl, Magn. Reson. Spectrosc. 18 (1986) 61-122 and in “Encyclopedia of Nuclear Magnetic Resonance”, Eds, D. M. Grant and R. K. Harris (Wiley, 1996) and Warren and Silver (Adv. Magn. Reson. 12 (1988) 247-384) describe that in addition to simple rectangular radio-frequency (rf) pulses with constant amplitudes and phases, composite and shaped pulses represent powerful tools for the manipulation of spins in NMR spectroscopy and imaging. In practice, both composite and shaped pulses are implemented as a sequence of rectangular pulses (with different amplitudes and phases). In the following, we will use the generic term “pulse” for both composite or shaped pulses. In the presence of experimental restrictions and errors such as maximum rf amplitudes and rf inhomogeneity, the attainable pulse performance of a pulse is limited by the pulse duration.
Depending on the task of a pulse in an experiment, there are desired and undesired single contributions of the system states. According to the present invention, desired single contributions of the system states are those single contributions of the system states that correspond to the system state demanded by the theory underlying a certain experiment. In contrast thereto, undesired system states are those single contributions of system states that differ from the system state demanded by the theory underlying a certain experiment.
Undesired single contributions to the system state comprise                phase errors        signals of nuclei having a certain resonance frequency (band selective pulses, solvent signal suppression)        decoupling sidebands        signal amplitudes deviating from a desired profile, e.g. wherein a desired profile can depend on frequency offset and/or the B1 field strength and/or coupling constants (comprising scalar coupling, dipolar coupling, quadrupolar coupling) and/or relaxation rates (comprising auto relaxation, cross relaxation, autocorrelated relaxation, cross-correlated relaxation) and/or diffusion rates and/or chemical exchange rates and/or chemical reaction rates and/or spatial coordinates, wherein a desired profile can be of any shape, e.g. constant signal amplitude, linearly dependent, quadratically dependent signal amplitude, polynomically dependent signal amplitude, exponentially or logarithmically dependent signal amplitude, a signal amplitude dependent on a random function,        signals of nuclei with undesired multiplicities        signals of undesired coherence orders (zero, single, double, triple quantum coherences and higher coherence orders)        signals from undesired coherence transfer pathways        signals having certain T1 and/or T2 relaxation timesor        signals of nuclei having direct coupling to a certain isotope (isotope filters).        
As an example, referring to a inversion pulse, when starting from equilibrium z magnetization, −z magnetization is a desired system state, all remaining components, for instance components along the +z axis as well as all transverse components are undesired system states in this example. Kobzar et al. (J. Magn. Reson. 170 (2004) 236-243 and J. Magn. Reson. 194, 58-66 (2008)) and Neves et al. (J. Magn. Reson. 181 (2006) 126-134) describe the determination of physical limits of pulse performance using methods of optimal control theory, as it has been described by Bryson and Ho (Applied Optimal Control, Hemisphere, Washington, D.C. (1975)). For example, for a given maximum rf amplitude and a desired bandwidth and robustness with respect to rf inhomogeneity, there exists a minimum pulse duration T* that is required to achieve a desired average fidelity or performance index. It is not possible for a single pulse to compensate its own imperfections to the desired degree if the pulse duration is shorter than T*.
Here, we show that pulse durations can be further reduced by allowing pulses to compensate each other's imperfections. In the following we will refer to this class of cooperatively acting pulses as COOP pulses. In multiscan experiments, for example, imperfections in individual scans are irrelevant if these imperfections cancel in the total accumulated signal. Keeler (Understanding NMR Spectroscopy, Wiley, Chichester, 2005), Bodenhausen et al. (J. Magn. Reson. 58 (1981) 370-388), Bain (J. Magn. Reson. 56 (1984) 418-427) Levitt et al. (J. Magn. Reson. 155 (2002) 300-306) describe phase cycling which is used in many multi-scan experiments for the suppression of artifacts or unwanted signals wherein in each scan, a sequence of identical (shaped) pulses is repeated, except for a systematic phase variation of the pulses (and the receiver). Here, we demonstrate that it is possible to improve the performance of pulse sequences by not only changing the overall phase of a given pulse in subsequent scans, but by cycling through a set of carefully designed COOP pulses which are in general not identical. Khaneja et al. (J. Magn. Reson. 172 (2005) 296-305) and Tosner et al. (J. Magn. Reson. 197 (2009) 120-131) describe the optimal-control-based gradient ascent pulse engineering (GRAPE) algorithm an adapted version of which can be used for optimizing highly compensating COOP cycles.
The optimization of a single (shaped or composite) pulse is usually conducted by the optimal-control-based gradient ascent algorithm GRAPE (“gradient ascent pulse engineering”).
Suppose for a given initial magnetization vector M(0) we are looking for a pulse of duration T that optimizes a defined performance index or quality factor φ, wherein we assume for simplicity (but without restricting the generality) that φ depends only on the final magnetization vector M(T). In the case of an excitation pulse, for example, we start with z magnetization, i.e. M(0)=(0, 0, 1)t. Skinner et al. describe (J. Magn. Reson 163 (2003) 8-15) how a simple quality factor can be defined as the x component of the final magnetization. A given pulse is fully characterized by the time-dependent x and y components νx(t)=−γ Brf,x(t)/2π and νy(t)=−γ Brf,y(t)/2π (or alternatively by the total if amplitude νrf(t)=√(νx2(t)+νy2(t)) and rf phase φ(t)=tan−1(νy(t)/νx(t)).
The pulse can be improved, if it is known, how the quality factor Φ reacts when the controls νx(t) and νy(t) are varied, i.e. if we know the gradients δΦ/δνx(t) and δΦ/δνy(t). These gradients can be approximated using finite differences.
Bryson et al. (Applied Optimal Control, Hemisphere, Wash. D.C. (1975)), Khaneja et al. (J. Magn. Reson. 172 (2005) 296-305), Tosner (J. Magn. Reson. 197 (2009) 120-134), Conolly et al. (IEEE Trans. Med. Imag. MI-5 (1986) 106115) and Skinner et al. (J. Magn. Reson. 163 (2003) 8-1) describe how the same high-dimensional gradients δΦ/δνx(t) and δΦ/δνy(t) can efficiently be calculated to first order based on principles of optimal control theory. This approach requires the calculation of the trajectory of the magnetization vector M(t), and of the so-called costate vector λ(t), for 0≦t≦T. Skinner et al. (J. Magn. Reson. 163 (2003) 8-15; J. Magn. Reson. 167 (2004) 68-74; J. Magn. Reson. 172 (2005) 17-23) describe that the desired gradients are approximated well by the x and y components of the cross product M(t)×λ(t):
                                          δΦ                                          δ                ⁢                                                                  ⁢                                                      v                    x                                    ⁡                                      (                    t                    )                                                              ⁢                                                                            =                                                                      M                  y                                ⁡                                  (                  t                  )                                            ⁢                                                λ                  z                                ⁡                                  (                  t                  )                                                      -                                                            M                  z                                ⁡                                  (                  t                  )                                            ⁢                                                λ                  y                                ⁡                                  (                  t                  )                                                                    ,                            (        1        )                                          δΦ                                    δ              ⁢                                                          ⁢                                                v                  y                                ⁡                                  (                  t                  )                                                      ⁢                                                                =                                                            M                z                            ⁡                              (                t                )                                      ⁢                                          λ                x                            ⁡                              (                t                )                                              -                                                    M                x                            ⁡                              (                t                )                                      ⁢                                                            λ                  z                                ⁡                                  (                  t                  )                                            .                                                          (        2        )            
For a spin with offset νoff, the effective field vector νe(t) is defined asνe(t)=(νx(t),νy(t),νoff)t,  (3)
Starting from the initial magnetization vector M(0)=Mi, the trajectory of the magnetization vector M(t) can be calculated by solving the Bloch equationsM(t)=2πνe(t)×M(t).  (4)
Rourke (Conc. Magn. Reson. 14 (2002) 112-129) describes modifications of the Bloch equations for taking into account relaxation, or radiation damping.
Here, for simplicity we assume that relaxation effects can be neglected, however if necessary they can be taken into account in a straightforward way. Khaneja et al. (J. Magn. Reson. 172 (2005) 296-305) and Gershenzon et al. (J. Magn. Reson. 188 (2007) 330-336) describe how relaxation effects can be taken into account in a straightforward way.
Skinner et al. describe (J. Magn. Reson. 163 (2003) 8-15) that if the pulse performance Φ depends only on the magnetization vector M (T) at the end of the pulse, the costate vector λ(T) is given by ∂Φ/∂M (T), i.e. the three components of the costate vector λ(T)=(λx(T), λy (T), λz(T))T are
                                                        λ              x                        ⁡                          (              T              )                                =                                    ∂              Φ                                      ∂                                                M                  x                                ⁡                                  (                  T                  )                                                                    ,                                            λ              y                        ⁡                          (              T              )                                =                                    ∂              Φ                                      ∂                                                M                  y                                ⁡                                  (                  T                  )                                                                    ,                                            λ              z                        ⁡                          (              T              )                                =                                                    ∂                Φ                                            ∂                                                      M                    z                                    ⁡                                      (                    T                    )                                                                        .                                              (        5        )            
For example, if the quality factor is simply the projection of the final magnetization vector onto a desired target state F, i.e.Φa=Mx(T)Fx+My(T)Fy+Mz(T)Fz,  (6)then the final costate vector is simply λ(T)=F. On the other hand, the quality to reach a target state F is defined asΦb=1−a1(Mx(T)−Fx)2−a2(My(T)−Fy)2−a3(Mz(T)−Fz)2,  (7)as has been described by Skinner et al. (J. Magn. Reson. 172 (2005) 17-23), then the resulting final costate vector is given by λ(T)=−(2a1(Mx−Fx), 2a2(My−Fy), 2a3(Mz−Fz))t. Here, a1, a2 and a3 represent the relative weights given to the desired match of the x, y, and z components of the magnetization vector and the target state. Skinner et al. (J, Magn. Reson. 163 (2003) 8-15; J. Magn. Reson. 167 (2004) 68-74; J. Magn. Reson. 172 (2005) 17-23;) and Gershenzon et al. (J. Magn. Reson. 188 (2007) 330-336) describe that the equation of motion for the costate vector has the same form as the Bloch equations (cf. Eq. (4)), i.e.{dot over (λ)}(t)=2πνe(t)×λ(t),  (8)and by propagating λ(T) backward in time we obtain λ(t) for 0≦t≦T.
Khaneja et al. (J. Magn. Reson. 172 (2005) 296-305) and Skinner et al (J. Magn. Reson. 163 (2003) 8-13) describe how robustness with respect to offset and rf inhomogeneity can be achieved by averaging the gradients over all offsets νoff and rf scaling factors s of interest. According to the present invention, averaging is to be understood as the calculation of the average value of individual system states, wherein the system states can be considered in the calculation of the average value with different weighting factors which can have both a positive as well as a negative sign and according to the present invention, the average value can be calculated according to the following methods comprising arithmetic, geometric, harmonic, quadratic or cubic averaging:
            arithmetic      ⁢                          ⁢      average      ⁢                          ⁢      value      ⁢              :            ⁢                          ⁢                        x          _                arithm              =                  1        N            ⁢                        ∑                      j            =            1                    N                ⁢                  x          j                                geometric      ⁢                          ⁢      average      ⁢                          ⁢      value      ⁢              :            ⁢                          ⁢                        x          _                geometr              =                            ∑                      j            =            1                    N                ⁢                  x          j                    N                  harmonic      ⁢                          ⁢      average      ⁢                          ⁢      value      ⁢              :            ⁢                          ⁢                        x          _                harmon              =          N                        ∑                      j            =            1                    N                ⁢                  (                      x            j                          -              1                                )                                quadratic      ⁢                          ⁢      average      ⁢                          ⁢      value      ⁢              :            ⁢                          ⁢                        x          _                quadr              =                            1          N                ⁢                              ∑                          j              =              1                        N                    ⁢                      (                          x              j              2                        )                              2                  cubic      ⁢                          ⁢      average      ⁢                          ⁢      value      ⁢              :            ⁢                          ⁢                        x          _                cubic              =                            1          N                ⁢                              ∑                          j              =              1                        N                    ⁢                      (                          x              j              3                        )                              3      wherein xj with jε{1, 2, . . . , N} corresponds to the elements that are averaged.
Starting from an initial pulse with rf amplitudes νx(t) and νy(t), the pulse performance can be optimized by following this averaged gradient. In the simplest approach, the gradient information can be used in steepest ascent algorithms, but faster convergence can often be found using conjugate gradient or efficient quasi-Newton methods that are also based on the gradients δΦ/δνx(t) and δΦ/δνy (t).
The method described above is limited to the optimization of one single pulse at a time.
The technical problem underlying the present invention is to provide an optimization method for being able to simultaneously optimize a group of pulses so that these pulses act in a cooperative manner.